Monte Carlo Methods: Lab 2

Point to note: an electronic copy of the Frenkel and Smit book is available through the library. This lab is based on case study 1 in chapter 3.4 of that book.

Lennard-Jones fluids

When computing the interactions between lots of bodies (atoms, molecules, planets, etc) we can either use the true potential or force between them, or we can approximate it with some potential (or force) that is easier (and usually cheaper) to calculate. The parameters of the potential can then be set to approximate the "real" features we're interested in.

In computational chemistry, one such approximation is the Lennard-Jones potential. Given two bodies separated by a distance $r$, the potential generated by those two bodies is

\begin{equation} U(r) = 4 \varepsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^{6} \right]. \end{equation}

Here $\varepsilon$ and $\sigma$ are parameters. When there are more than two bodies the total potential is the sum over all pairwise potentials.

In principle this generates a potential between particles that are separated by huge distances. Instead it is typical to truncate the potential: to pick a cut-off distance so that any particles separated by more than that distance do not contribute, and to correct for those small contributions.

Here we use a Lennard-Jones potential inside a box size $[0,L]^3$ with a cut-off $r_c = L/2$, with parameters set so that

\begin{equation} U = \begin{cases} 4 \left[ \frac{1}{r^{12}} - \frac{1}{r^6} \right] & r < r_c \\ 0 & r > r_c. \end{cases} \end{equation}

Include tail corrections (that is, additional energy and pressure terms resulting from the particles outside the cutoff radius) as

\begin{align} U^{\text{tail}} & = \frac{8 \pi \rho}{3} \left[ \frac{1}{3} \frac{1}{r_c^9} - \frac{1}{r_c^3} \right] \\ p^{\text{tail}} & = \frac{16 \pi \rho^2}{3} \left[ \frac{2}{3} \frac{1}{r_c^9} - \frac{1}{r_c^3} \right]. \end{align}

For each configuration we need to compute the pressure using

$$ \begin{equation} p = \frac{\rho}{\beta} + \frac{\text{Virial}}{V} \end{equation} $$

where

$$ \begin{equation} \text{Virial} = \sum_i \sum_{j > i} \vec{f}( \vec{r}_{ij} ) \cdot \vec{r}_{ij} \end{equation} $$

where, as usual, $\vec{r}_{ij}$ is the separation between the atoms, $\vec{r}_{ij} = \vec{r}_i - \vec{r}_j$, and the intermolecular force $\vec{f}$ is given by

$$ \begin{align} \vec{f}(\vec{r}_{ij}) &= - \nabla U \\ & = \begin{cases} 24 \left[ 2 \frac{1}{r^{14}} - \frac{1}{r^8} \right] \vec{r}_{ij} & r < r_c \\ \vec{0} & r > r_c \end{cases} \end{align} $$

Note that in the reduced coordinates $\beta = T^{-1}$.

Monte Carlo code

We will be using an $NTV$ approach, keeping the number of particles fixed ($N = 100$), the temperature fixed ($T=2$) and the volume fixed (indirectly, via the density $\rho = N / V = N L^{-3}$; use $\rho = a/10$ for $a = 1, \dots, 9$, but start by just considering the $a=1, 2$ cases). You will need to take at least $10,000$ steps for the larger values of $a$; $20,000$ is better, but in all cases you should test with a smaller number of particles and steps ($1,000$ may be sufficient for small values of $a$).

For reference we note the solutions, taken from Johnson, Zollweg and Gubbins for the pressures at $T=2$ are:

Efficiency

Note that the sum over all particles scales as $n^2$ where $n$ is the number of particles. As the number of steps the algorithm will need to take will also scale as $n$, this makes the number of calculations at least as bad as $n^3$. This is expensive; if you try the naive approach then you'll have difficulty using more than 50 particles in a moderate time.

Instead we can note that, at each stage, the algorithm will move only one particle. Therefore, if we store not just the locations of the particles but also their pairwise separations, at each step we will only have to modify a small number of the separations. So we can store $r^2_{ij} = \vec{r}_{ij} \cdot \vec{r}_{ij}$ only, for $j > i$, and when perturbing particle $k$ we only need to update the separations $r^2_{ik}$ for $i<k$ and $r^2_{kj}$ for $k<j$.

This should significantly reduce the number of calculations done in each step.

In addition, note that for reasonable behaviour the acceptance rate should be $\sim 40\%$. This depends on the fractional perturbation distance $\Delta$; values $\sim 0.4$ are reasonable when $\rho \sim 0.1$, but values $\sim 0.02$ are reasonable when $\rho \sim 0.9$.

Results

• Check that the energy has converged to a "constant" state.
• Plot a histogram of the energies to show that they follow the Boltzmann distribution.